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Abstract 

Model-based studies of auditory nerve responses to electrical stimulation can provide insight into the func- 
tioning of cochlear implants. Ideally, these studies can identify limitations in sound processing strategies and 
lead to improved methods for providing sound information to cochlear implant users. To accomplish this, models 
must accurately describe auditory nerve spiking while avoiding excessive complexity that would preclude large- 
scale simulations of populations of auditory nerve fibers and obscure insight into the mechanisms that influence 
neural encoding of sound information. In this spirit, we develop a point process model of the auditory nerve 
that provides a compact and accurate description of neural responses to electric stimulation. Inspired by the 
framework of generalized linear models, the proposed model consists of a cascade of linear and nonlinear stages. 
We show how each of these stages can be associated with biophysical mechanisms and related to models of neu- 
ronal dynamics. Moreover, we derive a semi-analytical procedure that uniquely determines each parameter in 
the model on the basis of fundamental statistics from recordings of single fiber responses to electric stimulation, 
including threshold, relative spread, jitter, and chronaxie. The model also accounts for refractory and summa- 
tion effects that influence the responses of auditory nerve fibers to high pulse rate stimulation. Throughout, we 
compare model predictions to published physiological data and explain differences in auditory nerve responses 
to high and low pulse rate stimulation. We close by performing an ideal observer analysis of simulated spike 
trains in response to sinusoidally amplitude modulated stimuli and find that carrier pulse rate does not affect 
modulation detection thresholds. 



1 Introduction 

Cochlear implants restore a sense of hearing to individu- 
als with severe to profound hearing loss by electrically 
stimulating the primary sensory neurons in the audi- 
tory pathway. A complete understanding of the trans- 
formation from electrical stimulation to neural responses 
would aid the design of improved cochlear implant sound 
processing strategies. Computational and mathematical 
models, especially those that are carefully constrained 
by available neurophysiological data, can play an essen- 
tial role in exploring this transformation. In addition to 
providing an efficient platform for testing and evaluat- 
ing stimulation paradigms, models provide quantitative 
tools for studying how neural mechanisms influence the 
transmission of sound information in the auditory sys- 
tem. 

In this paper, we develop a point process modeling 
framework that can be used to simulate the response 
of the auditory nerve (AN) to cochlear implant stim- 
ulation. The dynamical and stochastic features of the 
model are matched to statistics that characterize neural 
responses to single and paired pulses of electrical cur- 
rent - stimuli that are commonly used to characterize 
the response properties of AN fibers in animal models of 
electric hearing ( Hartmann et al.]|1984[ van den Honert 
Dynes||1996| |Bruce et al.||1999c 



and Mark||1992 


Litvak et al.|2003b Trevino et al.||2010 


Plourde et al. 


2011 


). Our approach is largely moti- 



vated by a specific family of point processes: gener- 
alized linear models (GLMs) ( |Paninski 2004 Truccolo 
et al.|2005[ ). GLMs have emerged as an essential tool for 
modeling physiological data and investigating the cod- 



ing and computational properties of neurons (Paninski 



2004 Paninski et al 



(Trevino et al. 2010 



20071, including auditory neurons 
Plourde et al. 2011). Moreover, 



they hold particular promise for applications to sensory 
neural prostheses because they have useful mathemati- 
cal properties that permit efficient parameter fitting to 



and Stypulkowski 



Miller et aJ.|T999] 



2000 



Miller et al. 



1984 



spike train data (Paninski 20041, they can be used to 
optimally decode spike trains ( Paninski et al.|2007 ), and 
they can be used in connection with real-time optimiza- 
tion methods to identify stimulation patters that control 
the timing of evoked spikes ( Ahmadi et al.pOlT ) . 

We apply our point process model to investigate dif- 
ferences between AN responses to low and high pulse 
rate electric stimulation. Interest in high pulse rate stim- 
uli arises from the intuitive notion that higher stimula- 
tion rates should provide greater temporal information 
to cochlear implant listeners and thereby improve speech 
reception. Indeed, computational modeling and neu- 
rophysiological studies have indicated that stimulating 
neurons with 5000 pulses per second (pps) pulse trains 



can desynchronize neural responses (Rubinstein et al. 



Shepherd and Javel|1999]|Cartee et al.| |1999||Litvak et al.|2001||2003a|. This may return AN ac 



2001a|b Cartee et al.|2006[ ). We go on tivity to a state that more closely resembles normal spon- 



to show that this model can provide insight into the re- 
sponses of neurons to extended pulse trains of electrical 
stimulation. Point process models also have mathemat- 
ical properties that facilitate analyses of auditory signal 
detection (Heinz et al. 2001a|b ), and we close by relat- 
ing results from our model to the psychophysical test of 
amplitude modulation detection. 

Our approach connects two prior modeling frame- 
works. The first is the stochastic threshold crossing 
model (Bruce et al. 1999a|c ). This model is a useful 
phcnomcnological representation of AN spiking. A num- 
ber of cochlear implant psychophysics experiments have 



taneous activity in the healthy cochlea ( Rubinstein et al 



1999 Litvak et al. 2003a). In related modeling and ex- 
perimental studies, high pulse rate stimulation has been 
shown to produce neural firing rates that more faith- 
fully reproduce the envelopc-s of temporally-modulated 



stimuli ( 


Litvak et al. 


2003c|b Chen and Zhang 


2007 


Mino 


2007). Psychophysical data from tests of ampli- 



tude modulation detection, however, have not yielded 
clear indications that high pulse rate stimulation im- 
proves the ability of listeners to detect temporal mod- 



been studied with this model ( 


Bruce et al. 


L999b 


Xu 




and Collins 


2007 Goldwyn et al. 


2010 


). Unfortunately, 



ulation ( 


Galvin III and Fu 


2005 


Galvin III and Fu 


2009 


). More g 



as stated by its creators, it lacks sufficient dynamical 
details to provide valid predictions of AN responses to 
high rate stimulation ( Bruce et al.||1999a |. Studying re- 



More generally, speech tests 
have reported mixed results as to whether higher pulse 
rate stimulation improves speech reception in cochlear 



sponses to high pulse rate stimulation is necessary to 
characterize contemporary cochlear implant stimulation 
strategies, and for reasons that we discuss in more detail 
below, is a primary focus of this work. We therefore ex- 
tend the stochastic threshold crossing model by turning 
to a more general class of point process models. 

Point processes describe spiking via an instantaneous 



implant listeners ( 


Brill et al. 


1997 


Kiefer et al. 


2000 


Loizou et al.||2000 


Vandali et al. 


2000 


Holden et al. 


2002 Friesen et al. 


2005). Computational modeling can 



firing rate that varies over time (|Perkel et al.|1967 


John- 


son 


199G 


Truccolo et al. 


2005). They have been fre- 



qucntly applied to model auditory nerve firing (Miller 



help trace the connections between the perceptual bene- 
fits, if any, of high pulse rate stimulation and the neural 
responses that it evokes. 

2 Method 

In this section we introduce the point process framework 
and discuss how it can be parameterized using commonly 
reported statistics of responses to single pulse and paired 
pulse stimuli. These include threshold, relative spread, 



1 



jitter, chronaxie, summation threshold, and two refrac- 
tory effects. We begin by introducing the response statis- 



tics that we used to construct the model (Sec. 2.1), fol- 
lowed by a general discussion of how these statistics can 



be extracted from a point process model (Sec. 2.2 1. We 



then introduce the specific model framework that we will 



use throughout the study (Sec. 2.3) and explain the pro- 



cedure by which each parameter in the model can be 



uniquely associated with a response statistic (Sec. 2.4). 
2.1 Response statistics 

A first-order description of neuronal excitability in re- 
sponse to electrical stimulation is the firing efficiency 
curve. This is an input/output function that relates the 
current level of a single pulse of current to the probability 
that the stimulus evokes a spike. As shown by Verveen 



in his pioneering recordings of green frog axons ( Verveen 



1960), the firing efficiency curve can take the form of a 
sigmoid function and be approximated by the integral of 
a Gaussian distribution. The median of the underlying 
Gaussian distribution represents the stimulus level that 
elicits a spike with probability one-half. This level is re- 
ferred to as the threshold of the neuron, which we denote 
by 9. A measure of the variability of spike initiation is 
the relative spread (RS), defined as the standard devia- 
tion of the underlying Gaussian distribution divided by 
its mean. Firing efficiency curves are frequently used 
to characterize the response of AN fibers to cochlear 
implant stimulation in electrophysiological experiments 



pynes|1996||Bruce et al.|1999c||Miller et al.|1999||Shep 



herd and Javel|1999 Miller et al.|2001b I and to parame 



et al. 


2002| Macherey et al. 


2007 


Imennov and Rubin- 


stein 


2009 


0' Gorman et al. 


2009 


. We will follow this 



approach and use the firing efficiency curve, in particu- 
lar the parameters 9 and RS, as the starting point for 
constructing our point process model. 

The firing efficiency curve summarizes the excitabil- 
ity of a neuron in response to a single pulse of stimu- 
lation. For longer pulse durations, the threshold cur- 
rent level is typically smaller, due to the capacity of the 
neural membrane to integrate charge over time. This 
dependence of threshold on pulse duration can be sum- 
marized by the chronaxie; the pulse duration at which 
the threshold level is twice what it would be for a much 
longer pulse duration (how long, and therefore reported 
values of chronaxie, vary in different experiments). This 
basic feature of membrane dynamics is absent from ear- 



lier stochastic threshold crossing models (|Bruce et al. 
1999a| |Litvak et al.|[2003b| 



The firing efficiency curve encapsulates variability in 
the initiation of spikes (or lack thereof), but an addi- 
tional source of randomness is in the timing of spikes. 
The standard deviation of spike times evoked by a sin- 
gle pulse is known as jitter. It is another statistic that 



has been widely studied in physiological and modeling 
studies of electrical stimulation of the auditory nerve 



van den Honert and Stypulkowski 



1984 



Miller et al. 



1999||Javel and Shepfaerd|2000| |Miller et al.|2001b||Car- 



teeet al.|2006[ ). Jitter can depend on pulse duration and 



pulse level, but in our model we will only use the value of 
jitter that is most commonly reported in the literature: 
the value measured for the same pulse duration and level 
that defines the threshold. 

In order for the model to generate realistic responses 
to high pulse rate stimulation, it must also include spike 
history effects. The first, and most essential, is a refrac- 
tory effect that reduces the excitability of the model neu- 
ron in the immediate aftermath of an evoked spike. This 
effect will be implemented as a transient increase in the 
threshold value 9 following a spike, similar to the method 
of Bruce and colleagues ( Bruce et al.||1999a ). The point 
process will also include a second spike history effect - a 



transient increase in RS immediately after a spike ( Miller 



et al. 2001a I. Based on simulation results from a bio- 
physically detailed computational model of an AN fiber, 
it has been hypothesized that this phenomenon is a sig- 
nature of the random activity of ion channels in a cell 



membrane known as channel noise (Matsuoka et al. 2001 



Mino and Rubinstein 20061. Modeling channel noise in 



the auditory nerve in the context of cochlear implant 
stimulation is important for a number of reasons. First, 
neurons in the deafened cochlea do not receive the typi- 
cal synaptic input from inner hair cells, and therefore ion 
channels may generate a dominant noise source in this 
stage of auditory processing. Second, channel noise may 
be a mechanism by which high pulse rate stimulation 
can desynchronize neural responses, thereby leading to 



improved encoding of electrical stimuli ( Rubinstein et al. 
19^9) |Litvak et al.||2003c|b| |Mino||2007[ ) 



Finally, we incorporate an additional feature into the 
model that is relevant for high carrier pulse rates - the 
capacity of a neuron to integrate consecutive subthresh- 
old pulses. In other words, the model will account for 
the phenomenon by which multiple pulses are more likely 
to evoke a spike than would be expected if each pulse 
acted independently on the neuron. This effect, known 
as summation or facilitation, has been observed in re- 
sponses to pairs of pulses with short interpulse intervals 
(Cartee et al. 2000 2006[) and to high rate pulse train 
stimuli (Heffer et al. 2010). As we will see, this form 



of interpulse interaction creates fundamental differences 
between responses to low and high rate pulse train stim- 
ulation. 

2.2 Point process theory 

In this section, we introduce a number of basic ideas from 
the theory of point processes and illustrate their connec- 
tions to the response statistics discussed above. A point 
process model is completely defined by its conditional 



2 



intensity function (Daley and Vere- Jones 2003), which 



can be interpreted as the instantaneous probability that 
a neuron will spike (Truccolo et al. 2005). We denote 



the conditional intensity function as X(t\I, H), where / 
is the stimulus applied to the neuron and H represents 
the past history of the neuron, both of which should be 
viewed as functions of time. A related and important 
quantity is the integrated intensity function: 



A{h,h\I,H) 



\{s\I,H)ds. 



From A(ti,t 2 \I,H), one can define the probability that 
there will be a spike in a time interval [ti,t 2 ,]- This 
function is known as the lifetime distribution function 
and is given by ( |Daley and Vere-Jones|[2003l ): 



L(h,t 2 \I,H) = l- e - A (*i^|/,ff)_ 



(1) 



The lifetime distribution function is central to our 
analysis because it forms the mathematical link between 
the point process model and the statistics that describe 
the responses of neurons to single and paired pulses. 
Consider, for instance, responses to a single pulse of cur- 
rent level / and some duration D. If we let t 2 in Eq. [I] 
become sufficiently large, then the lifetime distribution 
function defines the probability with which a single pulse 
of a certain current level / will ever evoke a spike. In 
other words, when viewed as function of /, it is equiva- 
lent to the firing efficiency curve. If we now change our 
perspective, and view Eq.[l]as a function of t 2 , for a fixed 
value of /, then this function gives the probability that 
a spike has been evoked before time t 2 . In this context, 
the lifetime distribution function describes the temporal 
dispersion of spike times, and can thus be used to cal- 
culate jitter. We now present how the single pulse and 
paired pulse response statistics discussed in Sec. |2.1| can 
formally be obtained using point process theory. 

We begin with a set of measures that are obtained 
from responses to single pulses, and we therefore assume 
that they are not influenced by spike history effects. This 
allows us to omit the term H for now and let / refer the 
current level of a single stimulating pulse. As introduced 
above, the threshold value is the current level at which 
the firing efficiency curve has a value of 1/2. From the 
definition of the lifetime distribution, must satisfy 



1 



= L(O,oc|0). 



Applying Eq. [T] this can be rewritten in terms of the 
integrated intensity function: 



A(O,oo|0) = log 2 



(2) 



To define RS, we generalize the traditional defini- 
tion that is based on assuming the firing efficiency curve 
is shaped like the integral of a Gaussian distribution 
( Verveen 1960). From Eq.[l] it is apparent that the firing 



efficiency curve for a point process model is not neces- 
sarily an integrated Gaussian function. It is, however, 
the cumulative distribution of some probability density 
function. To be concrete, we define an associated prob- 
ability density function as the derivative of the lifetime 
distribution function with respect to the current level /: 

/ J (0,oo|J) = ^:[A(0,oo|7)]e- A ( ' oo l 7 ) 

By analogy with the traditional definition, we let RS be 
the standard deviation of this associated density function 
normalized by its mean. 

To define chronaxie in terms of a point process model, 
we use Eq. [2] and compare threshold current levels for 
varying pulse durations. Let 9{D oa ) denote the threshold 
for a monophasic pulse of a long pulse of duration Doc . 
Then, the chronaxie is the phase duration D c for which 
the single pulse threshold 9(D C ) is twice the value of 
0(1)00). Specifically, D c must satisfy the following two 
conditions: 

9(D C ) = 29{D OQ ) 

A(0, oo|0(D c )) = A(0, oolflpoo)) = log(2) 

The final single pulse statistic we consider is jitter. 
As mentioned above, jitter can be obtained from the life- 
time distribution function when it is viewed a function of 
time. Following the standard practice, we will consider 
jitter in spike times evoked by a single pulse presented 
at the threshold level 0. Then the probability of at least 
one spike occurring in the interval [0,4], conditioned on 
the event that the stimulus produces a spike at any time, 
is twice the lifetime distribution function: 2L(0, 4|0). We 
define the associated probability density function by tak- 
ing a derivative with respect to time of Eq. [T] and see 
that jitter is the standard deviation of the following den- 
sity function for the probability that a spike occurs at 
time t: 



k(O,t\0) = 2A(O,4|0)e~ A ^ t l e >. 



(3) 



We now turn to responses evoked by pairs of pulses 
Using a summation pulse paradigm ( Cartee et al.||2000 



2006 ), we can quantify how threshold current levels change 



for a stimulus consisting of two pulses of equal current 
levels, separated by an interpulse interval (IPI) of vary- 
ing duration. Note that in experimental studies, the 
question asked is whether the pulse pair evokes a spike; 
the timing of the spike is not relevant. Assuming that 
the response to the pulse pair is not influenced by any 
prior spiking activity, the summation threshold satisfies 
the same threshold condition given by Eq. [2] the only 
changes are that the stimulus I is a pair of pulses and 
is the current level at which this pair of pulses evokes at 
least one spike on half of all trials, on average. 

Lastly, we observe how refractory effects can be in- 
corporated by allowing the single pulse threshold and 



3 



RS to depend on the time since the last pulse. In the re- as follows: 



fractory pulse paradi 


gm ( 


Cartee et al.|2000 


Miller et al. 


2001a 


Cartee et al. 


2006 


), a strong first pulse forces 



the neuron to spike and then a firing efficiency curve 
is measured from the responses to a second pulse pre- 
sented some time later. The mathematical relationships 
between that firing efficiency curve and the lifetime dis- 
tribution function remain unchanged, but they must be 
formally modified by adding the spike history term H. 

2.3 Model formulation 

We have presented the connection between the point pro- 
cess model and response statistics in a general manner. 
In this section, we make these connections more explicit 
by proposing a specific structure for the conditional in- 
tensity function. The model consists of a cascade of 
linear and nonlinear stages followed by a probabilistic 
spike generator. This structure is inspired by the pop- 
ular generalized linear model (GLM) class of point pro- 
cess neuron models. Fig. [I] illustrates the model fea- 
tures. The model differs from standard GLMs in several 
ways. These include a nonlinear stage that depends on 
the time since the last spike, an asymmetry in how the 
positive and negative phases of charge-balanced pulses 
are filtered, and a secondary filter that adds variability 
to simulated spike times. 

The action of the model, depicted schematically in 
Fig. [T] can be summarized as follows: an incoming pulse 
train I(t) is passed through stimulus filters K + (t) and 
K~(t), the outputs of the filters are recombined and used 
as input to a nonlinear function /(•). To incorporate re- 
fractory effects, the stimulus filters and the nonlincarity 
are all modified depending on the time since the last 
spike. An additional filter J(t) is included to add vari- 
ability to simulated spike times. The result of this chain 
of events is an instantaneous, history-dependent value for 
the conditional intensity function that defines the point 
process model: 

X(t\I, H) = [J* f{K+ *I+ + K~* r)](t), (4) 

where / + and J - represent the positive and negative 
portions of the input and * represents the convolution 
operator. 

A variety of methods exist for fitting GLMs to neural 
data including reverse correlation to white noise stimuli 



(Simoncelli et al. 20041 and maximum likelihood meth- 
ods (Paninski 2004). We pursue a different route and 
take advantage of the special structure of the proposed 
model and the mathematical relationships in Sec. |2.2| to 
develop a semi-analytical procedure that uniquely iden- 
tifies parameters in the point process model with spe- 
cific response statistics. For reasons of mathematical 
tractability and biological relevance, which we articulate 
more fully below, we define the components of the model 



f(v) = 
K+(t) = 

K(-t) = 

J(t) = 





if v > 





else 







-t/r K 


if t> 
else 


tk 




~t/T K 


if t> 
else 





-t/rj 


if t > 




else 



(5a) 
(5b) 
(5c) 
(5d) 



This model has five parameters: a, k, tk, ft, and tj. We 
next show how they are uniquely determined by the five 
response statistics discussed in Sec. |2.1| relative spread, 
threshold, chronaxie, jitter, and summation pulse thresh- 
old. We will also illustrate how the model can account 
for refractory effects by allowing the parameters a and 
K to depend on the time since the previous spike. 

2.4 Model parameterization 

We begin by discussing the response measures that do 
not depend on spike history, and thus neglect H for now. 
To simplify our presentation, we introduce some addi- 
tional notation. Let w(t) be the waveform of the stim- 
ulus; for instance, the waveform of a monophasic pulse 
is a rectangular step that reaches a maximum value of 
1. Alternatively, the waveform function for a biphasic 
pulse consists of two rectangular phases, one reaches a 
value of +1 and the other —1. We next define the filtered 
waveform function W(t), which represents the action of 
K + and K~ on w(t): 

W(t)= [ [K + (s)w + (t-s)+K-{s)w-(t-s)]ds, 
Jo 

where w + and w~ are the positive and negative parts of 
the waveform function. 

For the case that all pulses have identical current 
level /, for instance in single pulse or summation pulse 
experiments, the intensity function for the point process 
model has the compact form: 

\(t\I) = ( K I) a [J* W a ](t). 

The parameter k and the pulse current level / both fac- 
tor out due to our choice of a power law nonlinearity. In 
addition, the jitter filter J(t) can be ignored when con- 
sidering the integrated intensity function over all time. 



This follows from Fubini's Theorem (Jones 2001 ) and the 
fact that J °° J{t)dt = 1: 

poo pOO 

J*W a ](t)dt= / W(t) a dt. 
o Jo 

If we denote this integral as 

W(t) a dt. EE W a , (6) 
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Pulse 



Stimulus Nonlinearity Jitter Filter Probabilistic Spike 

Filters Spiking 



l(t) 



V K+ V 




Refractory 

Figure 1: Schematic diagram of the point process model. Current input to the model, shown here as a train of 
biphasic pulses, is passed through a cascade of linear filters and a nonlinear function. This produces a value - the 
conditional intensity - that defines the instantaneous probability of a spike, which is then used to generate random 
sequences of spike times. Previous spikes provide feedback that modulates the stimulus filters and the nonlinearity. 



then the lifetime distribution function for the time inter- 
val [0, oo) can be written as: 



distribution: 



L(0,oo\I) = l-exp(-n a I a W a ) 



(7) 



As noted above, this lifetime distribution function is, 
in the language of neurophysiology, the firing efficiency 
curve or input /output function of the neuron in response 
to a single pulse of applied current. Next, we use this 
equation to parameterize the model. 

2.4.1 Relative spread 

The choice of a power law nonlinearity significantly sim- 
plifies the parameterization process. The lifetime distri- 
bution function in Eq. [7] is a Weibull cumulative distri- 
bution function with shape parameter a and scale pa- 
rameter [k \/VVo] . We define the relative spread as 
the associated standard deviation divided by the mean. 



Using known expressions for these values we find ( Rinne 



2009): 



RS = 



T(l 



r 2 (i + i) 



-l, 



(8) 



where T(-) is the Gamma function. Remarkably, the re- 
lationship between RS and a does not depend on any 
other model or stimulus parameters. There is therefore 
a one-to-one relationship between the response statistic 
RS and the model parameter a. To fit the model, we 
invert this relationship to find a as a function of RS. In 
practice this must be done numerically. 

2.4.2 Threshold 

The threshold current level 6 is obtained by solving Eq.[5] 
which is equivalent to finding the median of the Weibull 



1 ./log 2 



(9) 



We invert this relationship in order to have an expression 
for the model parameter k as a function of the response 
statistic 8: 



1 

K = -r 



log 2 



(10) 



Note that k depends on the model parameters a as well 
as /?, and t k (through the definition of W„ in Eq. |6j. 
As we will see, the values of these parameters are inde- 
pendent of 9, and can therefore be viewed as (known) 
constants in this equation. 

2.4.3 Chronaxie 

In order to fit the model parameter t k , observe that the 
waveform function w(t) and thus the associated func- 
tion W a both depend on the pulse duration. We denote 
this dependence as W a {D). By definition, the chronaxie 
value is the pulse duration D c for which the threshold 
for a monophasic pulse is twice the value of the thresh- 
old in response to a monophasic pulse of a long duration 
Dqo (the exact value of varies across experimental 
studies). We can use Eq. [9] to define the threshold cur- 
rent level as a function of pulse duration, and obtain the 
relationship: 



1 I log 2 = 

K VW a (D c ) K 



2 „/ log 2 
W a (Ax>)" 



The factors of k cancel one another and the equation can 
be further simplified to: 



W a (D c ) 



= 2 C 



(11) 
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Conveniently, the solution to this equation depends only 
on W a and a. Now, a can be assumed to have a known 
value since it is uniquely determined by RS (see Eq. [8]). 
Furthermore, W a does not depend on (3 since chron- 
axie is defined with respect to monophasic pulses. As 
a consequence, the only undetermined model parame- 
ter in this equation is t k , which enters in the definition 
of W Q . There is no analytical way to identify t k as a 
function of chronaxie, but we can numerically evaluate 
yV a {Doo) and W a (D c ) and then use numerical root find- 
ing methods to determine the value of t k that satisfies 



the relationship in Eq. 11 for given values of chronaxie 
D c and reference pulse duration -Dqo. 

2.4.4 Summation threshold 

To parameterize the temporal summation properties of 
the model, we replicate the paired pulse experiments 



of Cartee et al. (20001 using pseudo-monophasic pulses 
( Cartee et al.|2000| " Cartee 2000). In these experiments, 
two charge-balanced pulses are presented at varying in- 
terpulse intervals (tin). The threshold current level for 
the combined response to both pulses, which we denote 
02 , is compared to the threshold current level for a single 
pulse presented in isolation, which we denote 6\. |Car-| 



tee et al. (2000) measured thresholds at three interpulsc 



intervals (100, 200 and 300 [is) and summarized the re- 
lationship between 9\ and 82 using the function: 



1 



1 



exp(-t/ P //T sum ). 



(12) 



The parameter r sum controls the length of time over 
which the neuron can effectively sum consecutive sub- 
threshold pulses. 

To translate this into the language of the point pro- 
cess model, we use Eq.[2]to obtain a relationship between 
the single and paired pulse thresholds. Let W a ,i denote 
the filtered waveform for a single pulse and W a ,2(tipi) 
be the filtered waveform for a pair of pulses separated 
by an interpulse interval of length i/p/. Then the ratio 
of the thresholds can be written as: 



W a , 2 (t/P/) 



(13) 



The model parameters that enter into this equation 
are a, t k , and /3 (since the pulses in this paradigm have 
negative phases). We treat a and t k as constants in 
Eq. [13] since they do not depend on /3 and can be de- 
termined at prior steps in the fitting procedure. We de- 
termine /3 by minimizing the mean square error between 
the two representations of 8 2 j8\ on the right hand sides 
of Eqs. [12] and [13] This must be done using numerical 
methods for integration and minimization. The parame- 
ter j3 will take a value between zero and one, with smaller 
values leading to greater subthreshold integration of con- 
secutive pulses. We are are not aware of experimental 



studies that that indicate how the temporal integration 
properties of the AN change immediately after a spike, 
so f3 is left as a constant that does not depend on spike 
activity. 

2.4.5 Jitter 

The remaining undetermined parameter is tj. To de- 
termine its value, we recall that jitter in the point pro- 
cess model is defined as the standard deviation of the 
probability density function in Eq. [3| Once again, we 
cannot obtain a simple analytical relationship between 
Tj and jitter, but we can use numerical integration and 
root finding methods to find the value of tj for which 
the standard deviation of the density function in Eq. [3] 
is equal to a measured value of jitter. In other words, we 
solve the minimization problem: 



argmm 



t 2 l t (0,t\8)dt- 



tl t (0,t\8)dt) -7 



where 7 is the experimentally observed jitter value. Al- 
though the density function l t (0, t\6) depends on all of 
the other parameters in the model, by this step in our 
parameterization method all other parameter values will 
have been specified and we can therefore treat them as 
(known) constants when determining the value of tj. 
Available data from cat AN fibers did not provide strong 
evidence that jitter changes significantly due to spike his- 
tory effects ( Miller et al.|200fa| ). Thus, the parameter tj 
does not depend on spike activity and is held constant 
throughout our simulations. 

2.4.6 Refractory effects 

Finally, we describe how spike history effects can be in- 
corporated in the model by allowing the parameters k 
and a to depend on the time since the last spike. The re- 



lationship between k and 9 in Eq. 10 shows that increas- 
ing n immediately following a spike allows the model 
to exhibit a common refractory property, whereby neu- 
rons are less excitable in the aftermath of a spike. Ad- 
ditionally, the relationship between a and RS in Eq. [8] 
shows that increasing a after a spike allows the model to 
exhibit an increase in RS during the refractory period. 
We choose the dynamics of these parameters to match 
the experimentally measured values reported in |Miller| 
et al. (2001a). This study used a paired pulse stimula- 



tion paradigm to probe the changes in threshold and RS 
due to spike time. Following an initial "masker" pulse 
that is sufficiently strong to always evoke a spike, there is 
an interpulse interval before a second pulse is presented. 
The current level of the second pulse is varied in order to 
measure the firing efficiency curve of the neuron within 
the refractory period. 



Miller et al. (2001a) defined the dependence of 8 and 



RS on the interpulse interval, which we denote as At, 



G 



with the following equations: 



4. Use a, t k , /3, and threshold value to determine the 



and 



9{At) 



RS{At) = 



1 _ e -(At-t 9 )/re 



RSn 



1 — e -(At-t RS )/T RS ' 



(14) 



(15) 



Here, 9q and RSq are baseline values obtained from sin- 
gle pulse responses. The parameter tg represents the 
absolute refractory period during which no spikes can be 
produced. We implement the absolute refractory period 
by holding the intensity function X(t\I,H) at zero un- 
til the elapsed time since the last spike exceeds tg. The 
time constant Tg quantifies how quickly the effects of a 
previous spike fade, and is therefore a measure of the 
relative refractory period. The parameters tas and trs 
play similar roles in defining the evolution of RS, as a 
function of the time since the last spike. 

To simulate the model with refractory effects, we let 
At in the above equations represent the time since the 
last spike. We then update the values of k and a at 
the onset of each pulse, based on the time since the last 
spike. These parameters are then kept at a constant 
value from the onset of the pulse until the onset of the 
next pulse in the stimulus, at which time their values are 
updated again to reflect the changed time since the last 
spike. This simplification, as opposed to allowing the pa- 
rameters to vary continuously, allows a straightforward 
parameterization of k and a using the same single pulse 
relationships with 6 and RS, respectively, which were 
derived above. The value of a must be obtained using 
numerical methods, but a formula for the value of k can 
be obtained by combining Eq. [10] and Eq. [14} 



«(Ai) = K (l-e- (At - te)/re ), 



(16) 



scale factor k in Eqs. 5b and 5c 



5. Use a, r K , j3, k and jitter value to determine the 
jitter time constant tj in Eq. |5d| 

6. Use the refractory function for RS in Eq. [15] and 
the relationship between RS and a to determine 
the dynamics of a after a spike. 

7. Use the refractory function for threshold in Eq. [14] 
and the relationship between threshold and k to 
determine the dynamics of k after a spike. 

Several comments are in order regarding the effect of 
pulse shape and duration on these parameters. Thresh- 



old should decrease with increasing pulse duration ( van d< in 



Honert and Stypulkowski 1984), this property is cap- 



tured for monophasic pulses by the chronaxie statistics 
and t k . However, thresholds also depend on pulse shape. 
In particular, thresholds measured with biphasic pulses 
are higher than those measured with monophasic pulse 
(IMiller et al.||2001b|). The model exhibits this prop- 



erty, but it is not quantitatively matched to experimental 
data. In contrast, the unique relationship between RS 
and a in Eq. [8] shows that, for this model, RS does not 
depend on pulse shape or pulse duration. This feature is 
consistent with data reported in ( Miller et al.||1999[ ), al- 
though others have suggested that RS may increase with 
pulse duration (Bruce et al. 1999b). The parameters (3 



where Ko is the baseline value that is determined from 
the single pulse threshold 6q. 

2.5 Summary of parameterization method 

The fitting sequence we have described above provides 
a semi-analytical method by which model parameters in 
Eqs. [5] can be uniquely determined based on single and 
paired pulse response statistics. This sequence can be 
summarized as follows: 

1. Use RS value to determine nonlinearity parameter 
a in Eq. |5a[ 

2. Use a and chronaxie value to determine the stim- 
ulus filter time constant t k in Eqs. |5b| and |5c| 

3. Use a, t k and the summation time constant to de- 
termine the value /3 in Eq. [5c] for scaling negative 
phases of the stimuli. 



and tj will both depend on the pulse shape and dura- 
tions in ways that we do not attempt to fit to physiolog- 
ical data. Ideally one would define the model with a set 
of data obtaining using a self-consistent set of stimula- 
tion parameters, and these same stimulation parameters 
would then be used to obtain additional physiological 
recordings or perform cochlear implant psychophysics 
experiments. We fit model parameters and compared 
subsequent simulation results to data from a variety of 
published sources; we will comment on variations in stim- 
ulation protocols and how they may impact our modeling 
results where appropriate. 

2.6 Numerical methods 

All simulations were performed using original computer 
code written in the Fortran programming language. Nu- 
merical quadrature to evaluate the stimulus and jitter 
filters was performed using the trapezoid method and 
a time step of 1 fis. Random numbers were generated 
using the Mersenne twister algorithm ( |Woloshyn]|1999[ ). 
The original Fortran code is available from the authors 
upon request and a more user-friendly version of the code 
written for Matlab (Mathworks, Inc.) can be down- 
loaded from the ModelDB website (accession number 



143760) or the website of one of the authors (http: / /www.cns.nyu. 
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K 12 



-k/2 



K I 



(refractory 



model neuron. We can reformulate the action of the 
stimulus filters K + and K~ in terms of a differential 
equation, and find that the state variable v(t) evolves 
according to: 



dv , _ . . . 



(17) 



P K I 



where the function g(I) is pictured in Fig.[2j\ and defined 
as: 



Current (I) 



9(1) 



n(At)I if I > 
Pk(M)I else 



B 




The function n(At), which depends on the time At since 
the last spike, is defined in Eq. [T6| 

The form of the differential equation for v(t) in Eq. 17 
reveals that its dynamics are similar to the widely-used 
leaky integrate and fire model ( Burkitt 2006 1 and related 



spike response models ( Gerstner and Kistlcr 



2002), a fact 



noted in previous presentations of GLMs ( Paninski|2004 



Paninski et al.|20"T0 1. The effect of spike history on g(I), 



Figure 2: Relationship between the point process model 
and a soft-threshold integrate and fire model. A: Illus- 
tration of the function g(I(t)) in Eq. 17 The parameter 
/3 controls the summation time constant of the model 
by decreasing the slope of g(I(t)) for negative currents. 
Refractory effects reduce the excitability of the neuron 
by decreasing the slope for all current values. B: Illus- 
tration of the probability density function for the noisy 
variable 77 in the stochastic threshold analogy (see text). 
Smaller values of a generate broader distributions, or 
equivalently a more variable spike initiation mechanism. 

3 Results 



which is determined by the dynamics of n(At) in Eq. [16} 
is to reduce the slope of this piecewise linear function, 
as shown in the transition from the blue to green lines 
in Fig. [2j\. This change in g{I) makes the model neuron 
less excitable immediately after a previous spike. This 
feature has also been included in spike response models 
( Gerstner and Kistlcr 2002 ) and integrate and fire mod- 



3.1 



els (Badel et al. 2008]). The asymmetry in g(I) enables 
the model to exhibit facilitation, or summation of con- 
secutive charge-balanced pulses. A modeling study using 
Hodgkin-Huxley type models has shown that summation 
of pseudomonophasic pulses, and in particular the time 
constant T sum in Eq. [l2j is determined by the dynamics 
of the m-gating variable ( Cartee|2000 ) . The asymmetry 
in g(I), therefore, can be viewed as a phenomenolgical 
approximation of the nonlinear process - mediated by 
Na + activation - by which the excitability of a neuron 
in response to consecutive charge-balanced pulses is in- 
creased beyond what would be expected if both pulses 
Connections to biophysical mechanism^ present ed in isolation from one another. 



and models 



In order to gain intuition into the structure and behav- 
ior of the point process model, we take a brief detour to 
explore its key features and parameters. To provide this 
intuition, we connect aspects of the point process model 
to biophysical mechanisms and more familiar mathemat- 
ical models of neurons. 

3.1.1 Subthreshold dynamics 

The result of applying the exponential filters K + (t) and 
K~(t) to the incoming stimulus I(t) can be equated with 
a dynamic state variable that we denote by v(t). This 
notation is meant to suggest, in a loose sense, that this 
state variable represents the membrane potential of the 



The parameter t k appears in Eq. [17] as the time scale 
of the dynamics of v(t). In the context of an integrate 
and fire model of a point neuron, this time scale is of- 
ten interpreted as the membrane time constant and re- 
lated to the passive integration properties of the neural 
membrane. In the context of cochlear implants, such a 
description overlooks the fact that intracochlear electric 
fields can interact with AN fibers at multiple, anatom- 
ically diverse regions of the neural membrane. For in- 
stance, membrane time constants of excitable nodes of 
Ranvier, somata, and unmyelinated portions of neurons 
may differ significantly (Rattay and Felix 2001). The 
parameter t k , therefore, must be interpreted as a sin- 
gle time scale that summarizes the combined integrative 
properties of a spatially extended AN fibers in an exter- 
nal electric field. 
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3.1.2 Spike generation 

In the second stage of the model, the power law nonlin- 



earity in Eq. 5a is applied to the state variable v(t) to 
generate the desired probability of spiking. In the GLM 
and related linear-nonlinear frameworks, this nonlinear- 
ity is typically static and often interpreted as a reduced 
description of the spike generating mechanisms of a neu- 
ron. We have introduced the innovation that the nonlin- 
earity depends on the time since the last spike in order to 
account for the observation that RS increases during the 
refractory period ( Miller et al. |[200ia I . Here, we explain 
the relationship between the nonlinearity and spike ini- 
tiation and show why a dynamic nonlinearity produces 
the desired change in RS. 

We begin by viewing the state variable v(t) along 
with the nonlinear function f(v) = v a as a definition 
of a point process conditional intensity function. The 
probability of observing a spike in a small time window 
St can then be obtained by approximating the lifetime 
distribution function in Eq. [T] We find that: 



P(spike; St) « 1 — e 



'St 



(18) 



We now consider an alternate point of view in which 
v(t) represents the subthreshold state of a neuron. We 
suppose that the probability that a spike occurs in the St 
time window is related to the distance between v(t) and 
some voltage threshold, which we denote 0. This idea 



is known as an escape noise model (Plesser and Gerst 



ner 



2000) and is closely related to stochastic threshold 
crossing models ( |Bruce et al.|1999b Litvak et al. 2003b). 



From this perspective, we can imagine that the noise free 
trajectory v(t) is modified by adding a random number rj 
that represents the noise in the spike generation process. 
The probability that a spike occurs in a St window can 
be written as: 



P(spike; St) « P{v + 7] > 0). 



(19) 



If we define the cumulative distribution function F. q for 
the noise variable n, then we can equate Eqs. [18] and [19] 
and get: 



F v (v) 



-(e- v ) a st 



By taking the derivative of this distribution function, we 
derive a probability density function for rj: 



f n {r l ) = a5t{Q-r 1 ) 



In Fig. [2)3, we illustrate the relationship between a 
and the probability density function of rj by plotting 
f v (ri) for several values of a. We computed these dis- 

The prob- 



_ ? /log2 



tributions using St — 1/is and = y —rr 
ability density function for rj becomes broad for small 
values of a. This shows how the parameter a character- 
izes variability in the spike initiation process. Moreover, 



this illustrates that allowing a to decrease immediately 
after a spike creates the desired behavior in the model - a 
spike generation mechanism that is more variable within 
the refractory period. To provide a biophysical interpre- 
tation for this parameter, we note that the random open 
and closing of ion channels in the cell membrane (known 
as channel noise) has been shown to determine the value 
of RS in simulation studies ( Rubinstein| 19 95 ) . Thus, we 
view rj as an abstract representation of channel noise. 

3.1.3 Spike time variability 

The final stage of the model is to apply the jitter filter 
J(i) to the output of the nonlinearity. As we will see 
below, this secondary filter is necessary in order for this 
point process model to have realistic amounts of spike 
time variability. J(t) is applied after the nonlinearity, 
which suggests that it represents a source of timing vari- 
ability subsequent to the spike generator in the cell. In 
the context of extracellular electrical stimulation of the 
AN, one possible source of this variability is spatial in- 
teraction among multiple possible spatial sites of spike 
initiation. For instance, in a simulation study of a spa- 
tially extended neuron with multiple nodes of Ranvier, 



Mino et al. (2004) showed that stimulating electrodes 



which are more distant from a model AN fiber produce 
more variable spike times. The more distant electric 
fields broaden the distribution of nodes of Ranvier at 
which spike initiation can occur, thus spatial interactions 
across the axon may represent one post-spike mechanism 
that adds jitter to spike times. 

3.2 Model parameterization 

With this understanding of the model features in hand, 
we proceed to parameterize the model following the steps 
outlined in Sec. |2.5| We defined parameters in the model 
using published data from single pulse and paired pulse 
recordings in cat. Table [l] summarizes the data values 
we used, their sources in the literature, and the param- 
eter values that we obtained from the fitting procedure. 
Whenever possible, we used mean values obtained from 
experiments that used biphasic pulses: these include the 



mean threshold, RS, and jitter values reported in (Miller 
et al.|2001b ). For chronaxie, we used the mean value re- 
ported by van den Honert and Stypulkowski ( 1984 ) for 
intracochlear electrical stimulation of cat AN fibers. The 
longest phase duration tested was 2000 /is. For the sum- 
mation time constant {r sum in Eq. 12), we did not find 



a mean value reported in the published literature, so we 
assumed a value that fell within the range reported in 



Cartee et al. (2006) for AN responses to an electrode 



placed in the scala tympani. The parameter values that 
describe the refractory effects on threshold and RS were 
obtained from Miller et al. (2001a I. Specifically, we used 
tg = 332/is and tq = 411/Lts, which were the mean values 
reported in this study. Mean values were not reported 
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Table 1: Neural response statistics and corresponding model parameter values 



Response statistic Value 



Reference 



Model Parameter Value 



Threshold 
Relative Spread 
Chronaxie 
Jitter 

Summation Time 



0.852 mA 
4.87% 
276 us 
85.5 lis 
250 lis 



( |Miller et al.]|2001b ) 

( Miller et aL|2001b ) 

( van den Honert"and Stypulkowski||l984 ) 
( Miller et al.||2001b ~ 
( Cartee et al]|2006| 



K 

a 

Tj 

P 



9.342 
24.52 
325.4 [ms 
94.3 lis 
0.333 
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Figure 3: Relationships between response statistics and 
model parameters. Following the procedure summarized 
in Sec. |2.5| we obtain unique relationships between re- 
sponse statistics (x-axis) and model parameter values (y- 
axis). At each step, a parameter value is chosen to fit the 
response statistic (black x , see Table [l] for exact values), 
and then this parameter value can be used to obtain a 
unique relationship between the next response statistic 
and model parameter value in the fitting hierarchy. The 



approximation used to determine a (Eq. 20 ) is shown as 
a dashed red line in the first panel. 

for tus and trs, so we estimated that tfts = 199 lis and 
T RS = 423/xs based on plotted data (Figure 8 in Miller 



et al. (2001a)). These experiments used monophasic 



stimuli, but it is straightforward to modify our refrac- 
tory model should data obtained using biphasic pulses 
become available. 

Fig. [3] depicts the relationships between the response 
measures shown on the x-axes and the model parameters 
shown on the y-axes. Each panel illustrates the unique 
relationships that arose by following the parameteriza- 
tion sequence summarized in Sec. |2.5| For instance, in 
panel A we show the dependence of the nonlincarity pa- 
rameter a on RS. Once we had fixed the value of a based 
on the desired RS value, we then computed how the stim- 
ulus filter time constant t k depended on chronaxie. This 
relationship is shown in panel B. One- by-one, we ob- 
tained values for each parameter, the values of which 
are marked by the x symbol in Fig. [3] and reported in 
Table Q] 

In order to introduce spike history effects on RS, 
we must recompute a at the onset of every pulse. In 
principle, we could use a numerical root finding method 
to invert Eq. [8j but to avoid this and thereby increase 
the computational speed of our simulations, we observed 
that the relationship between a and RS could be approx- 
imated by the power law: 



RS 



-1.0587 



(20) 



The combination of this equation for a and Eq. [15] pro- 
vided a simple rule for evolving a according to the time 
since the last spike. The approximation, which is shown 
with the red dashed line in the first panel of Fig. [3j is 
suitably accurate. 

3.3 Model predictions: single pulse stim- 
uli 

3.3.1 Firing efficiency curve 

We first simulated responses of the model to a single 
pulse of current. To be consistent with the experimental 
data from which we obtained values for threshold, RS, 
and jitter, we used a 40 lis per phase biphasic pulse stim- 
ulus as the input to the model. By varying the current 
level of the pulse, we could measure the complete firing 
efficiency curve, as shown in Fig. [4]A The probability of 
a spike is defined as the proportion of trials, out of a total 
of 5000 for each current level, on which the model gen- 
erated a spike. Simulation results are shown with black 
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Figure 4: Responses of the point process model to sin- 
gle pulse stimulation. A: Firing efficiency curves for the 
point process model (blue), an integrated Gaussian func- 
tion with the same mean and RS (red), and simulation 
results from the model (black x). B: Distributions of 
spike times obtained from 5000 simulations of the model 
with a jitter filter (blue) and without a jitter filter (red). 



x-marks. As expected, they line up with the analyti- 
cally derived firing efficiency curve obtained from Eq. [I] 
and plotted with a blue line. In order to compare this 
input/output function to the more standard functional 
form that is used in the stochastic threshold crossing 



model of Bruce et al. ( 1999c ), we show in red the integral 



of a Gaussian distribution with mean and standard de- 
viation consistent with the threshold and relative spread 
values that we used to fit the point process model. At 
the highest and lowest current levels, the firing efficiency 
curve for the point process model is slightly greater, but 
overall both methods produce similar spiking probabili- 
ties in response to a single pulse. 
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Figure 5: Responses of the point process model to the 
summation pulse paradigm. A: Simulated spiking prob- 
abilities (symbols) for three interpulse intervals: 200 (is 
(green), 500 (is (red), 1000 (is (blue). Solid lines are ana- 
lytically obtained from the lifetime distribution function 
in Eq. [I] Single pulse firing efficiency curve is shown 
in black for reference. B: Thresholds estimated from 
simulated firing efficiency curves and normalized by the 
threshold for two independent pulses are shown as x- 
marks for the three interpulse intervals. Facilitation oc- 
curs if this threshold ratio is less than one. Black curve 



is the analytically predicted result in Eq. 13 



of tj — 94.3/is, as determined by the parameterization 
method. This produces a distribution of spike times with 
standard deviation of 86 (is, consistent with the mean 



value in (Miller et al. 2001b) that was used to param 



eterize the model. In the absence of a second filtering 
stage, the model produces an extremely narrow distri- 
bution of spike times, as shown by the red distribution. 
The secondary filter is necessary, therefore, for the model 
to have realistic amounts of spike time jitter. 



3.3.2 Jitter 

Fig. [4j3 shows the distribution of 10000 spike times ob- 
tained from the point process model in response to a sin- 
gle biphasic pulse of phase duration 40 (is presented at 
the threshold stimulus level. The blue line shows shows 
results for the model with a jitter filter time constant 



3.4 Model predictions: paired pulse stim- 
uli 

3.4.1 Summation pulse paradigm 

To test how the neuron model responds to pairs of pulses, 
we first simulated the summation pulse procedure in 
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(Cartee et al. 2006). As inputs, we used two biphasic 



pulses of equal current level separated by an interpulse 
interval of either 200/xs, 500/zs, or 1ms. The probabili- 
ties of observing at least one spike in response to both 
pulses, for varying interpulse intervals and current levels, 
are shown as x -marks in Fig. [5]A_. These were computed 
from 5000 repeated simulations of the model. We also 
plotted the analytically obtained firing probabilities for 
pulse pairs, obtained from Eq. [I] as curves in panel A. 
The effect of summation is visible as a leftward shift of 
these curves at shorter interpulse intervals. As a ref- 
erence, the black line reproduces the single pulse firing 
efficiency curve from Fig. [4}\.. 

To summarize these results, we estimated summa- 
tion thresholds at each interpulse interval by fitting an 
integrated Gaussian function. We then normalized the 
estimated threshold values relative to the single pulse 
threshold and plotted the values in Fig. [5)3. The black 
curve represents the analytical prediction of the point 
process model, given by Eq. [13J for the threshold of two 
independent pulses. Facilitation occurs if the threshold 
ratio is less than one. In these simulations we see facilita- 
tion for interpulse intervals shorter than approximately 
1000 [is. The point process model, therefore, will exhibit 
facilitation for carrier pulse rates of roughly 1000 pps and 
above. 

3.4.2 Refractory pulse paradigm 

To probe the effects of spike history in the model, we fol- 



A 1 



lowed the experimental procedure in (Miller et al. 2001a 



and simulated responses to masker-probe pulse pairs. 
The current level of the first pulse was set to a very high 
value so that it always elicited a spike. The level of the 
second pulse was then varied in order to measure firing 
efficiency curves. We used pairs of biphasic pulses in or- 
der to be consistent with our previous simulations. Spike 
probabilities were defined as the proportion of trials (out 
of 5000 total) in which the second pulse produced a spike. 
Firing efficiency curves obtained from responses to the 
second pulse of current are shown in Fig.[6]A for three dif- 
ferent interpulse intervals (667, 1000, and 1500 fj,s). As 
the interpulse interval (and equivalently the time since 
the last spike) decreases, the spike history effect has a 
greater influence on the probability that the second pulse 
will evoke a spike. 

As discussed in Sec. |2.4.6[ there are two types of re- 
fractory effects in the model. The first is the standard 
refractory phenomenon whereby the model neuron is less 
excitable immediately after a spike. This is apparent 
in the rightward shift in the firing efficiency curves for 
shorter interpulse intervals in Fig. [6]A_. The second effect 
is that RS increases at shorter interpulse intervals. This 
leads to shallower slopes in the firing efficiency curves 
at shorter interpulse intervals, but the decibel scale in 
panel A obscures the change. There are some discrep- 
ancies between the simulation results (x -marks) and the 




B 



Q. 

© 
"D 



2 
o 



Current (dB re Single Pulse Threshold) 



1.5 



500 



1000 1500 

Time since spike ((is) 



2000 



CO 
N 



I 1-5 

Co 

£ 



500 1000 1500 

Time since spike (|xs) 



2000 



Figure 6: Responses of the point process model to the 
refractory pulse paradigm. A: Simulated spiking prob- 
abilities (symbols) for three interpulse intervals: 667 /is 
(green), 1000 fxs (red), 1500 (blue). Solid lines are an- 
alytically obtained firing efficiency curves obtained from 
lifetime distribution function in Eq. [TJ Single pulse 
firing efficiency curve is shown in black for reference. 
B: Thresholds estimated from the simulated firing effi- 
ciency curves and normalized by the single pulse thresh- 
old are shown as symbols for the three interpulse in- 
tervals. Black curve is Eq. [14] C: Relative spread es- 
timated from the simulated firing efficiency curves and 
normalized by the single pulse relative spread are shown 
as symbols for the three interpulse intervals. Black curve 
is Eq. pi 



predicted values from point process theory at the small- 
est interpulse intervals (green line). The differences arise 
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from the fact that the simulated spike times in response 
to the first pulse have a small amount of random vari- 
ability, whereas the analytical result is computed using 
the assumption that the time of the first spike is locked 
to the time of the onset of the first pulse. 

To summarize these effects, we estimated threshold 
and RS from the simulated firing efficiency curves by 
fitting integrated Gaussian functions. These values are 
normalized with respect to single pulse threshold and RS 
and shown in panels B and C, respectively. By construc- 
tion of the model, the values computed from simulations 
agree with the refractory functions Eq 
which we plot as black curves. We emphasize that a 
model with a static nonlinearity, namely one in which 
a does not depend on spike history, would produce RS 
values that decrease at shorter interpulse intervals. Al- 
ternatively, the stochastic threshold crossing model of 
assumes a constant RS and would 



The remaining panels of Fig. [7] show interspike inter- 
val distributions at the three pulse rates, where the cur- 
rent levels were set to evoke approximately 100 spikes 
per second. The distributions were obtained from his- 
tograms of 10000 spike times in time bins of 100 /is. The 
distributions obtained from responses to the 250 pps and 
1000 pps pulse trains, but not the high rate 5000 pps 
pulse train, show peaks that are clearly aligned to the 
interpulse intervals in the stimulus. This feature of the 
simulated interspike interval distributions is qualitatively 
similar to distributions recorded from cat AN fibers using 
14 and Eq. 15 the same pulse rates (Miller et al. 2008). For ease of com- 



3.5 



parison, we include examples of interspike interval data 
reported by Miller et al in the insets of each panel. The 
model is capable of reproducing some characteristics of 
the interspike interval distribution for this single neuron, 
although an important caveat is that the experimentally- 
measured interspike intervals were obtained under dif- 
ferent stimulus conditions and likely also different firing 
rates. We do not intend to claim that the model com- 

Model predictions: constant pulse train* lctel y reproduces all of the behavior of this or other 

cells. 



Bruce et al. ( 1999a 



produce a flat line in panel C. 



stimuli 



The results to this point have validated the fitting method 
and demonstrated that the model reproduces a range of 
measures that characterize responses to single pulses and 
pairs of pulses. Of greater relevance to cochlear implant 
speech processing strategies are the responses of neurons 
to extended trains of pulsatile stimuli. We first simulated 
responses to trains of biphasic pulses with phase dura- 
tion 40 /is and constant current levels, and sought to 
characterize how carrier pulse rate affects the sequence 
of evoked spikes. We relate our simulation results to 
relevant physiological data throughout. 

3.5.1 Firing rate and interspike interval distri- 
butions 

Fig. [7|^ shows how firing rates in the model increase 
with current level for three different pulse rates. The 
stimuli were one second in duration. Mean and stan- 
dard deviations (shown with error bars) were estimated 
from 100 repeated simulations. At the lowest pulse rate 
(250 pps, blue line), the interpulse interval is longer than 
the summation time scale as well as the relative refrac- 
tory periods for threshold and RS. Thus the firing rate 
curve follows directly from the single pulse firing effi- 
ciency curve in Fig. |4]A.. The firing rate saturates at one 
spike per pulse, in this case 250 spikes per second. The 
1000 pps and 5000 pps pulse trains have interpulse in- 
tervals of 1 ms and 200 /is, respectively, which are short 
enough for refractory and summation effects begin to im- 
pact the behavior of the model. Due to summation and 
the higher pulse rate, the same firing rate can be evoked 
by higher pulse rate stimuli using less current per pulse. 
Thus the firing rate curves are shifted leftward as the 
pulse rate increases. 



The results of these simulations suggest that spike 
time jitter, which was incorporated into the model on 
the basis of responses to single pulse stimuli, can account 
for distinctive features of the interspike interval distribu- 
tions and how they vary with the stimulus pulse rate. On 
the one hand, temporal variability in spike times is small 
relative to the interpulse intervals for the lower pulse rate 
stimuli, which leads to the periodic nature of these dis- 
tributions in Fig. [7p and D. On the other hand, this 
small temporal variability is sufficient to spread the sim- 
ulated spike times across the interpulse interval for the 
5000 pps, creating a more smoothly varying interspike 
interval distribution in panel B, consistent with the ap- 
pearance of the distribution shown in the inset. 

3.5.2 Synchronization of spike times 

Simulation and physiological studies have generated in- 
terest in the possibility that high rate pulse trains may 
desynchronize neural responses to cochlear implant stim- 
ulation ( |Rubinstein et al.|1999| jLitvak et al.|2003a[ ). We 
tested whether the model exhibited similar signatures of 
desynchronization by computing a synchronization in- 
dex, known as vector strength (IGoldberg and Brown 



1969), and then compared our results to physiological 
data reported in |Miller et al. ( 2008 1 . For a sequence 
of N spike times {ii}{ v , vector strength is defined with 
respect to a period T as: 



VS 
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.i=l 



cos(2nU/T) 
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.i=l 



sin(27rii/T) 



(21) 



VS takes values between and 1, with higher values be- 
ing interpreted as a more synchronized spike train. In 
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Figure 7: Firing rate and interspike interval distributions for three stimulus pulse rates. A: Firing rate as a 
function of current level per pulse. Stimuli were one second long trains of biphasic pulses with constant current 
level per pulse, the level of which is shown on the x-axis. Error bars indicate the standard deviation of the mean 
firing rate, estimated from 100 repeated simulations of the model. B-D: Interspike interval distributions estimated 
from 10000 interspike intervals with the pulse rates indicated in each panel. Current level per pulse was chosen so 
that the stimuli evoked approximately 100 spikes per second. Inset figures are single fiber recordings from cat AN 
fibers at the corresponding pulse rates reproduced from Fig. 1 
Springer Science+Business Media: J. Assoc. Res. Otolaryngol. 



Miller et al. (2008) with kind permission from 



these simulations, we use VS to measure the strength of 
phase locking to the period of the stimulus pulse train, 
and thus T represents the interpulse interval. In all sim- 
ulations, as well as the experimental data to which we 
compare our results, the stimulus is a train of biphasic 
pulses with a constant current level and a 40 /its phase 
duration. Simulation results for three pulse rates are 
shown in the left column of Fig. [HJ where each circle rep- 
resents the firing rate and VS values computed from the 
response to a 10 second long pulse train. Our results can 
be compared to VS values obtained from recordings of 
cat AN fibers responding to biphasic pulses trains pre- 
sented at the same pulse rates and pulse shape. These 

are 



(2008) 



data, which we reproduce from Miller et al. 
shown in the right column of Fig. [8] 

For the lowest pulse rate, VS values obtained from 
simulations exceeded 0.98 for all firing rates. This rep- 
resents near perfect phase locking to the 250 pps pulse 
train. VS systematically decreases with pulse rate - note 
the change of scale on the y-axes. Comparisons with the 



tween simulated and measured VS values. The primary 
discrepancy between the two sets of VS values is the 
large variability present in the data values, seen as scat- 
ter in the vertical direction of these plots. A likely cause 
of this difference is the fact that we simulated a single 
model neuron with a fixed set of model parameters. In 
contrast, the data were obtained from a sampling of 37 
AN fibers. Presumably each neuron had different in- 
trinsic properties that may even change further over the 
course of the experiment. Additionally, Miller et al sug- 
gested that some scatter in the VS values measured from 
responses to the 5000 pps stimulus may have been due to 
limitations in their ability to resolve the phase of spike 
times with respect to the 200 /is interpulse interval. 

On the basis of the lower VS values at higher pulse 
rates, one is tempted to conclude that these results re- 
veal the desynchronizing effects of high pulse rate stimu- 
lation. As noted in Miller et al. (2008), however, higher 



data of Miller et al. (2008) show good agreement be- 



VS values at higher carrier pulse rates may not indicate 
desynchronization since VS is computed with a different 
reference period for each stimulus, depending on the in- 
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Figure 8: Vector strength measured from responses to 
three pulse rates. Left column: VS computed from 
simulated responses to 10 second long pulse trains of 
biphasic pulses with equal current levels. Current levels 
were varied to explore a range of firing rates. Note the 
change of scale on the y-axis; VS values decrease as pulse 
rate increases. Right column: VS values reported by 



Miller et al. ( 2008 ) based on responses of 37 cat AN fibers 



stimulated with intracochlear electrodes using monopo- 
lar, biphasic pulse trains. Reproduced from Fig. 8 in 



(Miller et al. 2008) with kind permission from Springer 
Science+Business Media: J. Assoc. Res. Otolaryngol. 
The black diamonds do not represent median values, see 



Miller et al. (20081 for details 



terpulse interval. In particular, at higher carrier pulse 
rates the period T in Eq. [21] is smaller, and this can 
lead to lower VS values for spike trains with equivalent 
amounts of temporal dispersion. To test whether the 
low VS values computed in simulations using 5000 pps 
pulse trains do in fact indicate desynchronizing effects of 
high pulse rate stimulation, we examined the distribu- 
tion of spike times relative to the onset time of the pulse 
immediately preceding the spike. These spike time dis- 
tributions are shown in Fig. [9] and were estimated from 
10000 spike times using a 10/xs bin width. Panel A shows 
distributions for the 250 pps (blue) and 5000 pps (green) 
pulse trains for a current level set to evoke approximately 
100 spikes per second. In panel B, we show distribu- 
tions of spike times obtained from stimuli that caused 
the model to spike at 225 spikes per second. The spike 
time distributions are nearly identical at the lower firing 
rate indicating that the lower VS values obtained with 
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Figure 9: Further comparison of spike timing variabil- 
ity in response to low and high rate pulse trains. A: 
Distribution of times between simulated spikes and the 
onset of the previous pulse. Stimuli are trains of biphasic 
pulses of either 250 pps (blue) or 5000 pps (green) , with 
a constant current level per pulse. Current level was 
chosen so that firing rates of the model were approxi- 
mately 100 spikes per second for both stimuli. 10000 
simulated spike times were used to estimate the distri- 
butions, which were plotted using 10/xs bin window. B: 
Similar to A, with current levels increased so that firing 
rates in simulations were approximately 250 spikes per 
second. C: Similar relationship between jitter and firing 
efficiency in the point process model. Black curve rep- 
resents analytical calculation of jitter from lifetime dis- 
tribution function, blue circles mark the firing efficiency 
values, interpreted as average probability of a spike per 
pulse, corresponding to the pulse rates and firing rates 
in B and C. D: Relationship between and jitter and fir- 
ing efficiency in measurements of responses of cat AN 
fibers to monopolar, monophasic (cathodic) stimulation 
by an intracochlear electrode. Reprinted from Hearing 
Research ( |Miller et al.|1999[ ) with permission from Elsc- 



the 5000 pps stimulus do not reveal any desynchroniza- 
tion in this case. At the higher firing rate, however, the 
distribution of spike times in response to the 250 pps 
stimulus is considerably narrower than the distribution 
of spike times measured from the 5000 pps pulse train. 
One could interpret this as desynchronization, although 
it may be more accurate to say that the 5000 pps stimu- 
lus is maintaining a degree of spike time variability that 
is lost when the current level of the 250 pps stimulus 
is increased and the evoked firing rate approaches the 
maximal value of 250 spikes per second. 

We can explain the narrowing of the spike time dis- 
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tribution by considering the main source of spike time 
variability in the model. In Fig.[9p, we show the amount 
of jitter in simulated spike times, as a function of firing 
efficiency for single pulse stimulation. We obtained this 
curve directly from the density function for spike times 
in Eq.[3| For the 250 pps stimulus, firing rates of 100 and 
225 spike per second translate to firing efficiency values 
per pulse of 40% and 90%, respectively. The blue circles 
mark the location of these firing efficiency values. Jitter 
in the model decreases substantially at the higher firing 
efficiency value, which leads to the narrower spike time 
distribution in Fig. [9)3. In contrast, increasing the fir- 
ing rate from 100 to 225 spikes per second when using 
a 5000 pps stimulus translates to a small change in the 
firing efficiency, here interpreted as the average proba- 
bility of a spike per pulse, from 2% to 4.5% (green cir- 
cles). The model predicts, therefore, that high rate pulse 
trains can maintain temporal variability in spike times 
even at relatively high firing rates because the probabil- 
ity of a spike on any single pulse remains low. The key 
feature of the model that explains these results - the 
fact that jitter decreases with firing efficiency - has also 
been observed in the responses of cat AN fibers to elec- 



A Simulation 



trie stimulation using monophasic pulses (Miller et al. 
1999). These data are reproduced in Fig.[9t). They show 



a qualitative match with our point process predictions, 
and underline the conclusion that the desynchronization 
observed at high pulse rates largely follows from the use 
of weaker current impulses. 

3.5.3 Irregularity of firing responses 

An additional response statistic that has been used to 
characterize responses to high rate pulse trains is the 
Fano factor. This quantity measures irregularity in the 
number of observed spikes. It is defined as the variance 
in the number of spikes observed in a time window nor- 
malized by the mean value. A Fano factor of one, which 
is the value generated by any Poisson process model of 
spiking, is considered to signify a highly irregular spike 
generator and lower values indicate more regularity. To 
estimate Fano factor from the model, we simulated re- 
sponses to a one thousand second long train of biphasic 
pulses (40 /is pulse duration), using 250 pps and 5000 pps 
pulse trains. Ten estimates of Fano factor were obtained 
by subdividing these spike trains. Results are shown in 
Fig. [10]A., with error bars that represent standard error 
of the ten estimates of Fano factor for each stimulus con- 



dition. In Fig. 10 3, we reproduce a figure from (Miller 
et al.|2008[ ) in order to compare the behavior of the point 
process model to AN fibers in cat, also stimulated by 
biphasic pulse trains. These are median values obtained 
from recordings of 37 AN fibers in 8 deafened cats, as 
such they do not reveal the considerable variation around 
the medians present in the experimental data. 

Fano factors obtained from simulations in response 
to 250 pps pulse trains are shown by the blue curve in 
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Figure 10: Fano factor measured from responses to low 
(250 pps) and high (5000 pps) pulse rates. A: Simula- 
tion results using the point process model. Results for 
the low pulse rate stimulation (blue line) are predicted 
by the Fano factor for a binomial distribution (black line, 
see text for details). The gray line shows the correspond- 
ing prediction of a binomial distribution for the 5000 pps 
pulse train. The green line illustrates that simulated re- 
sponses to the 5000 pps pulse train can have larger Fano 
factors than responses to the 250 pps pulse train if the 
refractory period is shortened (re — 200 fxs). Fano factor 
values are estimated from 100 intervals of simulated spike 
trains, where each interval is 100 s long. Error bars rep- 
resent standard error in the mean of these estimates. B: 
Medians of Fano factors recorded from cat auditory nerve 



fibers, reproduced from Fig. 9 in (Miller et al. 2008) with 



kind permission from Springer Science+Business Media: 
J. Assoc. Res. Otolaryngol. 



Fig. \T0f i- and results for the high rate pulse trains are in 
red. The clear trend is a decrease in Fano factor as firing 
rate increases, consistent with the physiological data in 
panel B. As firing rates approach 250 spikes per second, 
the responses to the low rate stimulus saturate at one 
spike for every pulse, leading to a sharp decrease in the 
Fano factor, a trend also seen in the data. One sub- 
stantial discrepancy between simulation results and the 



median Fano factors reported in Miller et al. (2008) is 



that, in the data, Fano factors obtained from responses 
to high pulse rate stimuli were larger than those obtained 
from response to the low pulse rate stimuli for most fir- 
ing rates. The numbers above datapoints in Fig. [T0f3 
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indicate the multiplicative factor by which the high rate 
Fano factors exceed the low rate values. These experi- 
mental results supported the notion that high pulse rate 
stimulation can generate more irregular spiking activity. 

Can our model predict this phenomenon? As we show 
in greater detail in the Appendix, the effects of refractori- 
ness and interpulse interactions tend to decrease Fano 
factor; see also (Berry II and Meister 1998). When the 



response of a neuron to any pulse is independent of all 
previous stimuli and spike history, the activity can be de- 
scribed as a sequence of Bernoulli trials, where the prob- 
ability p that any one pulse evokes a spikes is equal to the 
average firing rate divided by the number of pulses. The 
total number of spikes, in this case, would be distributed 
binomially with mean Np and variance Np(l — p), where 
N is the total number of pulses. If we denote the aver- 
age firing rate by r and the pulse rate by p, then the 
Fano factor would have a simple dependence on p, and, 
consequently, the average firing rate r: 



Fi 



bino 



i-p = i 



(22) 



Our simulation results for the 250 pps pulse train are 
completely explained by this analogy to a binomial dis- 
tribution. Due to the relatively long interpulse interval 
of 4 ms, past spike and stimulus histories have no ef- 
fect on the model's response to subsequent pulses. The 
curve representing Ft,i nom i a i in Eq. [22] for the 250 pps 
pulse train is shown in black in Fig. |10| and follows the 
simulation results. At 5000 pps, there are strong effects 
of stimulus and spike history, so the binomial analogy is 
not expected to approximate the behavior of the model. 
However, it does provide an upper bound for the pos- 
sible Fano factors that the model can achieve, which is 
shown by the gray line. This illustrates that there is 
room for the model to generate higher Fano factor values 
in response to high pulse rate stimulation, but different 
parameter values must be chosen. To test the capacity 
for the model to produce higher Fano factors, we weak- 
ened the refractory effect by decreasing the time scale of 
the relative refractory period, Tg in Eq. [Mj to 200 ps. 
This value is still within the range observed in physio- 
logical recordings ( Miller et al.||2001b ). All other model 
parameters were kept the same. The results from sim- 
ulations of this model, shown in green, illustrate that it 
can achieve higher Fano factor values, although the dif- 
ference between Fano factors at the low and high pulse 
rates remains larger in the data. 

We were able to explore the relationship between re- 
fractory dynamics and Fano factor in the model more 
fully by developing a discrete-time Markov chain approx- 
imation to the point process model. The Markov chain 
framework extends the Bernoulli process analogy (intro- 
duced above for low pulse rate stimulation) for stimuli 
in which spike history and interpulse interaction effects 
are present. We used the Markov chain approximation 
to obtain estimates of firing rate and Fano factor and 



explore the full range of possible behaviors in the model. 
See the Appendix for further details. 

3.6 Model predictions: amplitude-modulated 
pulse train stimuli 

Modern cochlear implant speech processors provide tem- 
poral information to cochlear implant listeners by mod- 
ulating the current levels of pulses over time. It is im- 
portant, therefore, to explore how the model responds to 
non-constant pulse trains. In order to relate model re- 
sults to available physiological and psychoacoustic data, 
we simulated responses to sinusoidally amplitude-modulated 
pulse trains. Specifically, as inputs to the model we used 
trains of biphasic pulses with the current level of the n th 
pulses defined by the equation: 

In = ^(1 + msin(27rt„/ m )). 

/„ is the current level of the pulse with onset at time t n , 
and the parameters m and f m parameterize the modula- 
tion depth and modulation frequency, respectively. As in 
previous sections, these pulses were charge balanced and 
had phase duration of 40 /is. To compare the model with 



available physiological data (Litvak et al. 2003b) and a 
recent modeling study ( O'Gorman et al.||2010 ), we used 
a 5000 pps carrier pulse rate modulated at a frequency 
of 417 Hz. We simulated two mean current levels, one 
that evoked approximately 50 spikes per second when the 
pulse train was unmodulated and a second that evoked 
approximately 100 spikes per second. The duration of 
all pulse trains used in simulations was ten seconds, and 
the results from ten repeated simulations for each stim- 
ulus condition were used to compute standard errors in 
estimated response statistics. 

In Fig.[y}\, we show simulated firing rates, as a func- 
tion of modulation depth. Firing rates increase with 
modulation depth for m greater than approximately 5%. 
We compared these results to cat single fiber data mea- 



sured by Litvak and colleages ( Litvak et al. 2003b ) , which 
we have reproduced in Fig. 11 H. The simulation results 



qualitatively match the experimental data, especially when 
compared to the two fibers that have lower firing rates 
(circle and filled square in C) . The main discrepancy be- 
tween the simulated and recorded firing rates is that the 
model does not exhibit increased spiking until the modu- 
lation depth increases beyond 5%, whereas the recorded 
firing rates increase at modulation depths as low as 1%. 

To quantify the sensitivity of spike timing to the pe- 
riod of the modulated waveform, Litvak et al computed 
The relevant period in this case is the 

Simula- 



VS as in Eq. 21 



21 



period of modulation, so T=2.4ms in Eq. 
tion results using the point process model are shown in 
Fig. 11 3 and the experimental measurements of |Litvak| 



et al. (2003b I are reproduced in panel D. The model 



again qualitatively captures the increase in VS with in- 
creasing modulation depth, although with important quan- 
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Figure 11: Responses to sinusoidally amplitude- 
modulated pulse trains with 417 Hz modulation fre- 
quency. A,C: Simulation results obtained from re- 
sponses to ten second long biphasic pulse trains with 
40 lis pulse durations and modulation frequency given 
on the x-axis. Two mean current levels were used, one 
that produced a mean firing rate of 50 spikes per sec- 
ond for unmodulated input (blue) and a higher value 
that produced a mean firing rate of 100 spikes per sec- 
ond (red). Error bars represent standard error in the 
mean of ten repeated simulations. VS is computed with 
reference to the period of modulation. B,D: Firing rate 
and VS measures obtained from five cat AN fibers stimu- 
lated by intracochlear monopolar cochlear implant stim- 
ulation using biphasic pulse trains with duration of 25 
jj,s per phase. Reprinted with permission from Litvak| 
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et al. ( 2003b I , copyright Acoustical Society of America. 



titative differences. In particular, the experimental data 
show that AN fibers have exquisite sensitivity to weak 
modulations, producing VS values of approximately 0.5 
for modulation depths of only 0.5%. As shown by O'Gorman 

ing of pulses (see Fig 



Figure 12: Histograms of response to 250 pps (left col- 
umn) and 5000 pps (right column) pulse trains modu- 
lated at 75 Hz. Modulation depth increases from 1% 
(top row) to 10% (bottom row). Histograms were esti- 
mated from 10000 repeated simulations of responses to 
one cycle of trains of biphasic pulses, and plotted with 
50 fis time bins. Current levels were set separately for 
each pulse rate so that unmodulated pulse trains pre- 
sented would evoke approximately 50 spikes per second. 



et al. (20101, this sensitivity to weak modulations may 



be a consequence of nonlinear mechanisms in spike gen- 
erators that are accounted for in Fitzhugh-Nagumo dy- 
namical models, but appear to be lacking in this point 
process description. 

3.7 Application to amplitude modulation 
detection 

In this final set of simulations, we explore how the tem- 
poral information in spike patterns can be quantified in 
a way that enables simulation results to be interpreted 
in the context of psychoacoustic experiments studying 
the ability of cochlear implant listeners to detect modu- 



lation at varying carrier pulse rates ( 


Galvin III and Fu 


2005 


Pfingst et al. 


2007 


Galvin III and Fu 


2009 


). We 



begin with a simple observation - a 250 pps carrier, due 
to the tight phase-locking of spikes relative to the tim- 



, does not appear to evoke 
AN activity that represents a modulated envelope. This 
point is illustrated in the left column of Fig. [l2j Here, 
we show histograms of spike times obtained from one 
cycle of a pulse train modulated at 75 Hz. From top 
to bottom, the modulation depth is increased from 1% 
to 10%. The current level is fixed so that 50 spikes per 
second are evoked, on average, in response to unmodu- 
lated stimuli. Our intuition tells us that this sparse and 
punctate pattern of spikes will not effectively transmit 
information about a slowly-varying envelope. Nonethe- 
less, these responses still carry some information about 
modulation depth, as evidenced, for instance, by the in- 
creasing height of the second peak as modulation depth 
increases. 

In the right column of Fig. 12 we show responses to 
a cycle of a 5000 pps pulse train modulated at 75 Hz. 
These histograms indicate that the high rate carrier ap- 
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Figure 13: Percent discrimination of a sinusoidally 
amplitude-modulated pulse train. Three discrimination 
measures are were used to discriminate unmodulated 
pulse train from a modulated pulse train on the basis of 
simulated spike trains (see text for details). Both stim- 
uli were one second long trains of biphasic pulses, AOfis 
per phase and current level was varied for each carrier 
pulse rate so that the unmodulated pulse train always 
evoked approximately 50 spikes per second, on average. 
Modulated pulse trains had 75 Hz modulation frequency 
and 1% modulation depth. Percent correct values were 
computed from 1000 repeated simulations of the model, 
error bars represent standard errors in the mean of ten 
repeated estimates of percent correct. Chance level is 
50% correct. 



pears to provide a more complete representation of the 
envelope of the modulated stimulus. In contrast to the 
low rate carrier, then, we may expect that this carrier 
would provide greater temporal information regarding 
the presence of a sinusoidal modulation. In other words, 
we may expect improved amplitude modulation detec- 
tion if an observer (for instance higher processing centers 
in the auditory pathway) had access to this pattern of 
spikes as opposed to those in the left column. 

To test these statements quantitatively, we used an 
ideal observer analysis to simulate modulation detection 
based on spike trains generated by the point process 
model. A general result of point process theory allows us 
to use the conditional intensity function X(t\I(t),H(i)) 
in Eq. [4] to compute the log-likelihood of observing a 
sequence of spike times {ti}f , conditioned on the input 
stimulus I {t) ( |Daley and Vere- Jones 2003 Paninski 2004 



Truccolo et aI1[2005| ): 

L({t i }^\I)=J2^g(X(t i \I(t i ),H(t i ))- 



N 



X(t\I(t),H(t))dt, (23) 



where s is the total duration of the stimulus. To an- 
alyze how temporal modulations are encoded in spike 
trains, we simulated two spike trains: one in response 
to an unmodulated pulse train and one in response to a 
modulated pulse train with 75 Hz modulation frequency 
and 1% modulation depth. Both stimuli were one second 
long and consisted of biphasic pulses that were 40/is per 



phase. After computing the likelihood function (Eq. 23 ) 
for all possible pairings of spike trains and stimuli, we 
used a likelihood ratio test to discriminate between the 



two spike trains ( 


Green and Swets 


1966 


Pillow et al. 


2005 


Goldwyn et al.|2010|. If the true pairings of stim- 



ulus and response produce the highest likelihood values, 
then the ideal ideal observer is said to correctly detect 
the modulated stimulus. By repeating this procedure, 
we estimated percent correct detection values. We also 
varied the carrier pulse rate from a low carrier pulse rate 
(250 pps) to a high carrier pulse rate (5000 pps) to in- 
vestigate our previous observation that low pulse rate 
stimulation produces an (apparently) incomplete repre- 
sentation of the modulated waveform. 

The results of this analysis are shown by the red line 
in Fig. \13\ where error bars indicate the standard error 
in the mean of 10 computed values of percent correct. 
The x-axis shows the carrier pulse rate used as the input 
to the point process model. At all pulse rates tested, for 
this controlled spike rate of 50 spikes per second, this 
maximum likelihood discrimination procedure produced 
correct discrimination on approximately 80% of trials. 
These simulations did not reveal any strong changes due 
to the carrier pulse rate. Despite the visible differences 
in Fig. |12| the spike patterns in response to modulated 
stimuli are equally detectable regardless of the carrier 
pulse rate. 

To further characterize the information in these spike 
patterns, we also computed discrimination measures us- 
ing decision rules that selected the modulated pulse train 
on the basis of which spike train response had a higher 
spike count (blue line) or a higher vector strength with 
reference to the 75Hz modulation (green line). From 
Fig. |11| we can see that spike count and vector strength 
increase with modulation depth, so simulated spike trains 
can encode the presence of modulation information us- 
ing either of these cues. For the small modulation depth 
used in these simulations, the spike count discrimination 
rule had near chance performance at just over 50% cor- 
rect. The vector strength discrimination measure pro- 
duced a higher percent correct, although it was still sub- 
stantially below the value obtained using the spike train 
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likelihood technique discussed above, indicating that vec- 
tor strength as a measure of phase locking does not com- 
pletely describe how temporal modulations are expressed 
in the precise sequence of spike times. Similar to the 
likelihood function discrimination measure, neither spike 
count nor vector strength discrimination predicted per- 
cent correct values that depended strongly on the carrier 
pule rate. 

To summarize these simulations of modulation detec- 
tion: we found that modulation detection based on three 
different decision rules did not vary across a range of 
carrier pulse rates. This ideal observer analysis suggests 
that modulation detection does not require peripheral 
spike patterns that fully represent the sinusoidal enve- 
lope of these pulse trains. 

4 Discussion 

4.1 Review of main findings 

In this study, we have used point process theory to derive 
mathematical expressions for statistical measures of neu- 
ral excitability that are commonly used to characterize 
the response properties of AN fibers to electrical stimu- 
lation. Furthermore, we proposed an explicit model with 
features that could be related to biophysical mechanisms 
in a phenomenological sense. These included tempo- 
ral filtering that represented subthreshold dynamics of 
the membrane potential, a nonlinearity associated with 
spike generating processes, and a secondary filter that 
accounts for variability in spike timing. 

An essential feature of the proposed modeling frame- 
work is that all parameters can be determined from re- 
sponses to single and pairs of electric pulses. The model 
is minimal in the sense that each model parameter is 
uniquely identified with a single statistic reported in 
physiological experiments, as illustrated in Fig. [3] More- 
over, the point process model incorporated dynamical 
and stochastic properties that are relevant to high pulse 
rate stimulation. Specifically, the relative spread of the 
model depends on the time since the last spike in a man- 
ner consistent with data reported in ( Miller et al.|2001a ), 
and subthreshold pulses presented in rapid succession 
could combine to increase the excitability of the model 
neuron, a facilitation phenomenon observed in electrical 
stimulation of AN fibers ( |Cartee et al.|2000[ |2006| |Heffer 
etH||20l0| . 



The model produced interspike interval distributions that 



qualitatively agreed with data reported in (Miller et al. 
2008 ) (see Fig. [7|) , and also provided important insight 



into how temporal jitter in spike times could account 
for synchronization to pulse trains presented at multi- 
ple pulse rates (see Fig. |8| . Overall, it appears that the 
proposed point process framework provides a satisfactory 
description of AN spike trains, although there are impor- 
tant shortcomings and potential extensions to consider, 
which we discuss below. 

4.2 Relation to past modeling studies and 
future directions 

The point process model represents an extension to pre- 
vious simplified models of AN spiking. In particular, 
the model developed in Bruce et al. ( 1999a|c ) and re- 



lated stochastic threshold crossing models (Litvak et al 



2003b) by including temporal integration, facilitation, 
spike history dependent relative spread, and a realis- 
tic amount of spike time jitter. The simplest dynamical 
models that describe AN responses to electric stimula- 



tion are integrate and fire models (Stocks et al. 2002 



Chen and Zhang 2007). As explained in Sec. 3.1 the 



point process model can be interpreted as a modified in- 
tegrate and fire model with escape noise. One significant 
difference between our approach and standard integrate 
and fire models is that we have introduced an asymmetry 
in the otherwise linear subthreshold dynamics, enabling 
the model to exhibit facilitation in response to charge 
balanced biphasic pulses. An advantage of the point pro- 
cess model over standard integrate and fire models is that 
spiking probabilities can be computed relatively simply 
using the conditional intensity function and the corre- 
sponding lifetime distribution functio n in Eq. [Tj One 



does not need to estimate hazard rates ( Plesser and Ger- 
stner|2000 ) or solve difficult first passage time problems. 



We have used this mathematical tractability to derive 
a parameterization method on the basis of physiological 
data. 

We suggest that this feature is an important strength 
of our model. Cochlear implants stimulate thousands of 
AN fibers, and response statistics such as threshold, rel- 
ative spread, and jitter can vary widely across the pop- 



ulation of AN fibers, see for instance (Miller et al. 1999 



2001b). Constructing populations of AN fibers with a 



representative distribution of thresholds has been done 



The construction of the model ensured that it would 
reproduce the measures of threshold, relative spread, jit- 
ter, chronaxie, as well as changes in threshold and rela- 
tive in the context of response to pairs of pulses separated 
by an interpulse interval. To further test whether our 
model generalizes, we simulated responses to extended 
stimuli of greater relevance to cochlear implant stimu- 
lation: pulse trains with constant current per pulse and 
current levels that were sinusoidally amplitude-modulated. 



using a biophysically detailed model ( jlmennov and Ru- 
binstein 2009), but the point process framework pro- 



vides an explicit method for controlling a number of 
key response properties of model neurons. In addition 
to facilitating future studies that consider the response 
properties of populations of AN fibers, this paramet- 
ric control can also be used to investigate which fea- 
tures may impact sensory perception in a significant way. 
For instance, physiological data from rat AN fibers sug- 
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gest that long-term deafness can increase absolute refrac- 



tory periods and excitability thresholds (Shepherd et al. 



2004 ) . Such changes can be incorporated into the model 



via straightforward modifications of parameter values. 

Although the model accurately predicted responses 
to extended pulse trains in several ways, it was unable 
to quantitatively predict the precise dependence of Fano 



factor on pulse rate and firing rate (Fig. 10) and the 



exquisite sensitivity of AN fibers to weak modulations 
depths as small as 0.5% (Fig. [8]). An explanation for 



these results can be found in O' Gorman et al. (2010), 



which showed that the Fitzhugh-Nagumo exhibits a dy- 
namical instability when driven with 5000 pps stimuli, 
and that this mechanism can account for experimentally 
measured values of Fano factor and strong phase locking 



to weak modulations reported by Litvak et al. (2003b). 



Future work could seek to synthesize these developments 
by developing models that contain this dynamical mech- 
anism and can also be both parameterized to additional 
physiological data and directly used in likelihood-based 
discrimination studies. One possible approach for im- 
proving the sensitivity of the model to small modulation 
is to apply an additional filter to the stimulus that acts as 
a differentiator on the sinusoidal evelope. Filters could 
also be included, for instance, to account for resonator 
dynamics that may affect AN responses to cochlear im- 
plant stimulation (Macherey et al. 20071. The analogy 



of subthreshold dynamics with escape noise discussed in 



Sec. 3.1 can provide useful insight into the effects of ad- 
ditional filters on the point process model. 

An alternative approach, that would still be grounded 
in the point process framework, would be to fit a model 
directly to a richer set of stimuli that are more relevant to 
the study of cochlear implant speech processing strate- 
gies. We have pointed out that the form of the model 
presented here, a cascade of linear and nonlinear stages, 
is similar to standard GLMs (jPaninskil 120041) . GLMs 



have been shown to capture the activity of sensory neu- 
rons in retina with a high degree of temporal precision 
and can be fit from sets of spike trains recorded in re- 



sponse to arbitrary time- varying stimuli (Pillow et al. 
2005 ) . They have also been recently applied to auditory 



nerve recordings for acoustic stimulation (Trevino et al. 



2010 Plourde et al. 2011). An important direction for 



future work, therefore, would be to fit GLMs to single 
unit recordings of AN fibers using stimulus sets that are 
clinically relevant to cochlear implant speech processing. 
One could also pursue a model-based approach and fit 
GLMs to biophysically detailed models of AN responses 



to cochlear implant stimluation (Imennov and Rubin- 



stem 



2009 Woo et al. 2010). The resulting point pro- 



cess descriptions would not necessarily be constrained by 
data from single pulse and paired pulse stimuli, but the 
mathematical relationships in Sec. |2.2| could still be used 
to analyze the resulting models. 

Incorporating firing rate adaptation is another im- 



portant direction for future models of AN responses to 
cochlear implant stimulation. Recordings from AN fibers 
in deafened cats have shown that adaptation affects neu- 
ral responses to amplitude-modulated cochlear implant 
stimuli (Hu et al. 2010[) and differentially affects responses 



to low and high rate pulse trains (Zhang et al. 20071 



GLMs can exhibit firing rate adaptation through the ad- 
dition of longer lasting spike history effects. Compu- 
tational modeling can provide further insight into the 
origins and coding consequences of adaptation for stim- 
ulation strategies ( Woo et aL]|2010 |. 



4.3 Implications for high pulse rate stim- 
ulation strategies 

A motivation for this work was to develop a model that 
could accurately predict responses to high carrier pulse 
rate stimulation. We have done this by including dynam- 
ical and stochastic features that reflect how past stim- 
ulation and past spiking activity can influence the ex- 
citability of AN fibers. Our simulation results in Fig. 12 



past modeling studies (Rubinstein et al. 1999 Mino et al. 



2002 Chen and Zhang 2007 ) , and physiological evidence 
~( Litvak et al. 2003c|b I suggest that responses to high 
pulse rate stimulation may represent the envelope of tem- 
porally modulated stimuli with greater fidelity than re- 
sponses evoked by low pulse rate stimulation. Moreover, 
one would intuitively expect high pulse rate stimulation 
to provide greater temporal information to cochlear im- 
plant listeners and thereby improve speech perception 
and psychophysical performance on tests of temporal res- 
olution. The fact that, to date, no psychoacoustic studies 
have shown clear evidence that high pulse rate stimu- 
lation improves speech reception or amplitude modula- 
tion detection poses an intriguing challenge to our under- 
standing of the connection between AN spiking activity 
and auditory perception. 

To connect AN spiking activity to psychoacoustic ex- 
periments, one can leverage the mathematical theory of 
point processes which provides access to tools from signal 
detection and information theory (Heinz et al. 2001a|b 



Goldwyn et al. 2010). In Sec 



3.7| we used the likeli- 
to 
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hood function of the point process model (Eq 
simulate modulation detection for a 75 Hz sinusoidally 
amplitude-modulated pulse train. In human listeners, 
modulation detection at this frequency is correlated with 
speech perception ( Won et al.|2011 ), so it is of interest to 
understand how neural responses transmit the temporal 
information in these stimuli. 

The simulated spike trains evoked by low pulse rate 
stimulation were equally discriminable to those evoked 
by high pulse rate stimulation in our simulation of am- 
plitude modulation detection (Fig. 13). These results 



suggest the possibility that cochlear implant listeners 
could successfully discriminate modulated and unmod- 
ulated stimuli in the context of psychoacoustic exper- 
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imcnts, even when patterns of evoked neural activity 
in the auditory periphery provide distorted or incom- 
plete representation of the stimulus envelope. This ob- 
servation is consistent with psychoacoustic evidence that 
modulation detection in human listeners does not im- 



prove with high carrier pulse rate stimulation (Galvin 



HI and Fu||2005] |Pfingst et aL]|2007| |Galvin III and Fu 
2009). Moreover, it highlights an essential challenge for 



improving listening outcomes for cochlear implant users 
in order to identify the relative benefits (or shortcom- 
ings) of high pulse rate stimulation, research should seek 
to identify psychoacoustic measures that reveal the per- 
ceptual consequences of the distinct patterns of neural 
activity evoked by low and high rate pulsatile stimula- 
tion. 
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Appendix: Markov chain approxi- 
mation to the point process model 

Definition of the Markov chain for con- 
stant pulse train stimulation 



In Sec. |3.5.3| we discussed how the firing rate and Fano 
factor produced by low pulse rate, constant current level 
stimulation can be estimated for the point process model 
by making an analogy to a Bernoulli process. In this Ap- 
pendix, we generalize this type of approximation for the 
case of high pulse rate, constant current level stimula- 
tion by using a Markov chain model to account for the 
effects of past spikes and stimulus history. Results from 
Markov chain theory allow us to characterize the range 
of firing rates and firing irregularity (Fano factor) that 
the model can produce for this class of stimuli. 

As in the Bernoulli process analogy, we simplify the 
problem by associating to each pulse a probability that 
the neuron spikes in response. This probability is given 
by the lifetime distribution function in Eq. [TJ where the 
integral is evaluated over a duration of one interpulse in- 
terval. This probability value depends on the history of 
past pulses and past spikes. If we neglect precise spike 
timing, and presume that all spike times coincide with 
the onset of a pulse, then we can approximate the prob- 
abilities at hand via a sequence of values {p n (I)}, where 
the subscript n represents the number of pulses that have 
elapsed since the last occurrence of a spike. In the ab- 
sence of history effects, for instance for low pulse rate 
stimulation, p n (I) is identical for all pulses since the in- 
terpulse interval is long relative to summation and re- 
fractory effects. 



To account for history effects, we define a discrete 
time Markov chain. The states of this chain, which 
we denote s n , represent the number of pulses that have 
elapsed since the previous spike. There are two possible 
transitions away from each state. If a spike occurs, then 
the chain returns to s±. The transition probability of 
this event is p n (I)- If no spike occurs, then the state of 
the chain advances by one to s n +i. The probability of 
this event is 1 — p n (I). In practice, there is some number 
of pulses beyond which the refractory and summation 
effects no longer alter the probability of a spike, so we 
can limit the number of states in the chain to some fi- 
nite number TV. A schematic illustration of the resulting 
Markov chain is shown in Fig. [14]A An example of the 
sequence of transition probabilities {p n (I)} is shown in 
Fig. [I4j3 for a 5000 pps stimulus that produces a firing 
rate of 100 spikes per second. There are 25 total states in 
this chain (the x-axis) , this indicates that history effects 
in the model do not persist beyond ~ 5ms. We therefore 
used a 5 state Markov chain for 1000 pps stimuli and a 
2 state Markov chain for 250 pps stimulation. 

Calculation of the firing rate and Fano fac- 
tor 

Results from the theory of discrete-time, discrete-space 
Markov chains allow us to analytically compute the firing 
rate and Fano factor of the Markov chain approximation, 
for arbitrary pulse rates and current levels. To do this, 
we use the fact that the Markov chain has a stationary 
distribution n which gives the long time probability that 
the Markov chain will be in each state. The stationary 
distribution can be obtained by computing the eigenvec- 
tor associated with the (unique) eigenvalue that takes 
the value one ( |Kemeny and Snell|1960 ). In other words, 
the stationary distribution ir solves 

TtM = 7T, 

where M is the transition matrix for the Markov chain. 
The first element of tt, which we denote by 7r 1; represents 
the long-term proportion of time steps in which the chain 
is in state si. Since the chain only returns to state s\ if 
there has been a spike in response to the previous pulse, 
7Ti can be used to compute the firing rate. In particular, 
if we denote the pulse rate of the stimulus by p, then the 
firing rate (in spikes per second) is given by: 



(24) 



Fig. 14 



compares this Markov chain approximation for 
firing rate (black lines) to firing rates obtained from sim- 
ulations of the point process model (colored error bars) . 
The simulated firing rates are the same as shown in 
Fig. [7]A The fact that the Markov chain approxima- 
tion accurately predicts firing rates for the point process 
model at low pulse rates is expected because the spik- 
ing can be described by a Bernoulli process in this case. 
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Fi gure 14: Markov chain approximation to the point process 
model. A: Schematic of the Markov chain, with p n representing 
the probability that the neuron spikes n pulses after the previous 
spike. B: Example of p n as a function of s n for 5000 pps stim- 
ulation set to evoke 100 spikes per second (blue). Corresponding 
piecewise linear approximation (red), see text for details. C: Fir- 
ing rates of the point process model predicted by the Markov chain 
model. D: Fano factor of the point process model predicted by the 
Markov chain model. E: Dependence of the Fano factor of the 
piecewise linear approximation to the Markov chain model on the 
number of initial and recovery stages. Contours show Fano factor 
values with firing rate set to be 100 spikes per second and pulse 
rate of 5000 pps. Gray area indicates regions where Fano factor 
is higher at 5000 pps than the value obtained at 250 pps, and 
X marks the position of the point process model with parameter 
values given in Table [T] 



Importantly, however, the Markov chain also accurately 
approximates firing rates for the high pulse rate stimuli. 
This framework, therefore, adequately captures the in- 
fluence of spike history and interpulse effects in the point 
process model and can be used to estimate the firing rate 
of the point process model via Eq. [24] 

Next, we show how the Markov chain approximation 
to the point process model can be used to compute the 
Fano factor of the point process model for pulse train 
stimuli with constant current strength. If we let N(T) 
be the number of spikes produced over a period of time 
of length T, then the Fano factor is defined as: 



F 



Var[7V(T)] 



The denominator is given via the stationary distribu- 
tion, as discussed above. In particular, E[iV(T)] = ■n\n vi 
where n p is the total number of pulses in the stimulus. 
To estimate the variance of N(T), we first define the 
characteristic matrix ( Kemeny and Sneli]|1960 ) 



Z = [Id - M + Afoop 1 

where Id is the identity matrix, M is the transition ma- 
trix for the Markov chain, and is a matrix whose 
rows are the vector n. We then appeal to the law of 
large numbers for Markov chains, which states that the 
Umiting variance of iV(T)/ v /n^ is 7Ti(2Zn — tti — 1) (Ke- 
meny and Snell|1960 ), where Z\\ is the first entry of the 
characteristic matrix. 

To compute the Fano factor, we then take the ratio 
of the variance of N(T)/ \ffip to its expected value and 
find that the Fano factor is: 



F = 2Z n - tti - 1. 



(25) 



Fig. [I4j3 shows values of the Fano factor for 250 pps and 
5000 pps stimulation computed using the equation for 
a range of firing rates (black lines). Fano factor values 
computed for the point process model are shown with 
colored error bars, and are the same as those in Fig.[7}3. 
Overall, there is close agreement between the Markov 
chain approximation given by Eq. [25] and the Fano fac- 
tor values computed from simulations of the point pro- 
cess model. In sum, the Markov chain approximation 
accurately captures both the mean firing rate and the 
normalized variance of the spike count for a range of 
stimulus levels and pulse rates. 

Piecewise-linear approximation to the Markov 
chain 

As discussed in the text, physiological data suggest that 
Fano factors measured from responses to 5000 pps pulse 
trains may be higher than those measured from responses 
to 250 pps pulse trains ( Miller et aL||2008 1 . We were in- 
terested in identifying whether the point process model 
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could exhibit the same behavior, so we used the Markov 
chain to systematically explore the space of possible Fano 
factors that the model could produce. To do this, we 
first observed that for the 5000 pps pulse train, the p n {I) 
curves can be roughly caricatured as having three stages: 
an initial stage where the probability of spiking is near 
zero, a recovery stage where p n {I) increases in a rela- 
tively linear manner with n, and a final stage where the 
p n {I) is nearly constant with n. This piecewise linear 
approximation is shown in red in Fig. [14)3. We therefore 
characterized the functions p n (I) by two parameters, one 
that describes the number of initial stages and one that 
describes the number of recovery stages. We then swept 
across this space of two parameters to explore the range 
of possible Fano factor values for a fixed firing rate. 

An example of this procedure is shown in Fig. |14E , 
where firing rate is set to be 100 spikes per second and 
the pulse rate is 5000 pps. We show Fano factor values 
computed using the Markov chain model with the piece- 
wise linear p n (I), and indicate in gray the regions of 
parameter space in which Fano factors are higher when 
stimulating the neuron with 5000 pps than with 250 pps 
pulse trains. The black x indicates the location of the 
model with parameter values given in Table [I] In order 
to generate Fano factors that are larger for responses to 
5000 pps stimuli than 250 pps stimuli, either the num- 
ber of initial stages or the number of recovery stages 
must be decreased. One mechanism in the model that 
controls the number of initial and recovery stages is the 
time scale of the threshold refractory effect, Tg in Eq.fTi] 
Shortening the refractory period reduces the number of 
initial and recovery stages, and thus our analysis of the 
Markov chain model provides a theoretical basis for the 



observation in Fig. 10 that the point process model with 
a smaller Tg will have higher Fano factor for 5000 pps 
stimulation than 250 pps stimulation. The number of 
initial and recovery stages is also influenced by the time 
scale of stimulus integration, t k in the stimulus filters 
K + and K~ in Eqs. 5b and 5c The Markov chain ap- 



proximation is therefore a useful tool for connecting be- 
tween spike train statistics of the point process model to 
the biophysical properties (refractory effects, membrane 
time constant, e.g) that are represented by parameters 
in the model. 
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